# With Village Fixed Effects ---------------------------------------------------

# LPG Connection (Binary)
model.dd.vfe.uselpg.scst = fixest::feols(uselpg ~ scst*wave + age + gender +
                                           religion + education + hhsize +
                                           lnexp + bplaay | village, 
                               data = pool, weights = pool$weights)

# LPG Home Delivery (Binary)
model.dd.vfe.convenience.scst = fixest::feols(convenience ~ scst*wave + age + gender +
                                           religion + education + hhsize +
                                           lnexp + bplaay | village, 
                                         data = pool, weights = pool$weights)

# Grid Connection (Binary)
model.dd.vfe.usegrid.scst = fixest::feols(usegrid ~ scst*wave + age + gender +
                                            religion + education + hhsize +
                                            lnexp + bplaay | village, 
                                          data = pool, weights = pool$weights)

# With Household Fixed Effects -------------------------------------------------

# Regression Models for Energy Access Response Variables
# Robust standard errors are clustered at household level.

# LPG Connection (Binary)
model.dd.hfe.uselpg.scst = fixest::feols(uselpg ~ scst*wave + hhsize +
                                           lnexp + bplaay | hh, 
                                         data = pool, weights = pool$weights)

# LPG Home Delivery (Binary)
model.dd.hfe.convenience.scst = fixest::feols(convenience ~ scst*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool, weights = pool$weights)

# Grid Connection (Binary)
model.dd.hfe.usegrid.scst = fixest::feols(usegrid ~ scst*wave +  hhsize +
                                            lnexp + bplaay | hh, 
                                          data = pool, weights = pool$weights)

## Model and coefficient lists for village and household fixed effects combined ##
models.vfe.hfe.dd_access_scst <- list(model.dd.vfe.uselpg.scst, 
                                  model.dd.vfe.convenience.scst,
                                  model.dd.vfe.usegrid.scst, model.dd.hfe.uselpg.scst,
                                  model.dd.hfe.convenience.scst, 
                                  model.dd.hfe.usegrid.scst)

esttex(models.vfe.hfe.dd_access_scst, 
       se = "cluster", 
       digits = 3,
       dict = c(scst = "SC/ST", wave = "2018", uselpg = "LPG Connection",
                usegrid = "Grid Connection", convenience = "LPG Home Delivery",
                village = "Village", hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_access_scst.tex"))
       
# With village and household fixed effects and two interaction terms -----------

# LPG Connection (Binary)
model.dd.hfe.uselpg.sc_st = fixest::feols(uselpg ~ sc*wave + st*wave + hhsize +
                                           lnexp + bplaay | hh, 
                                         data = pool, weights = pool$weights)

# LPG Home Delivery (Binary)
model.dd.hfe.convenience.sc_st = fixest::feols(convenience ~ sc*wave + st*wave + hhsize +
                                                lnexp + bplaay | hh, 
                                              data = pool, weights = pool$weights)

# Grid Connection (Binary)
model.dd.hfe.usegrid.sc_st = fixest::feols(usegrid ~ sc*wave + st*wave +  hhsize +
                                            lnexp + bplaay | hh, 
                                          data = pool, weights = pool$weights)

## Model and coefficient lists for village and household fixed effects combined ##
models.hfe.dd_access_sc_st <- list(model.dd.hfe.uselpg.sc_st,
                                       model.dd.hfe.convenience.sc_st, model.dd.hfe.usegrid.sc_st)

esttex(models.hfe.dd_access_sc_st, 
       se = "cluster", 
       digits = 3,
       dict = c(sc = "SC", st = "ST", wave = "2018", uselpg = "LPG Connection",
                usegrid = "Grid Connection", convenience = "LPG Home Delivery",
                hh = "Household"),
       keep = c(":", "2018"),
       fitstat = c("r2"),
       fixef_sizes = T,
       float = F,
       replace = T,
       file = c("./Manuscript/Tables/dd_access_sc_st.tex"))

# Visualisation of linear combination of interaction terms ---------------------

# Replicates the stata lincom command, using robust clustered standard errors
# for our specific case of one interaction term (wave) and two binary 
# variables (sc and st) interacted separately

lincom_custom <- function (fit, interactionterm, var1, var2) {
  
  sumvar1interaction <- fit$coefficients[[interactionterm]] + fit$coefficients[[var1]]
  sumvar2interaction <- fit$coefficients[[interactionterm]] + fit$coefficients[[var2]]
  
  vcov <- summary(fit)$cov.scaled
  
  sevar1interaction <- sqrt(
    fit[["se"]][interactionterm]^2 +
      fit[["se"]][var1]^2 +
      2 * vcov[interactionterm, var1])
  
  sevar2interaction <- sqrt(
    fit[["se"]][interactionterm]^2 +
      fit[["se"]][var2]^2 +
      2 * vcov[interactionterm, var2])
  
  lincom_df <- tibble(.rows = 3,
                      "term" = c(interactionterm,var1,var2),
                      "coef" = c(fit$coefficients[[interactionterm]],
                                 sumvar1interaction,sumvar2interaction),
                      "se" = c(fit[["se"]][[interactionterm]],
                               sevar1interaction, sevar2interaction))
  
  return(lincom_df)
  
}

uselpg_lincom <- lincom_custom(model.dd.hfe.uselpg.sc_st, "wave", "sc:wave", "wave:st") %>% 
  mutate(model = "LPG Connection")

convenince_lincom <- lincom_custom(model.dd.hfe.convenience.sc_st, "wave", "sc:wave", "wave:st") %>% 
  mutate(model = "LPG Home Delivery")

usegrid_lincom <- lincom_custom(model.dd.hfe.usegrid.sc_st, "wave", "sc:wave", "wave:st") %>% 
  mutate(model = "Grid Connection")

lincom_df <- rbind(rbind(uselpg_lincom, convenince_lincom), usegrid_lincom) %>% 
  mutate(term = case_when(
    grepl(term, pattern = "sc") ~ "SC",
    grepl(term, pattern = "st") ~ "ST",
    TRUE ~ "All Others"
  ))

lincom_df %>% 
  ggplot(aes(model, coef, colour = term, shape = term,
             ymin = coef - 1.96 * se,
             ymax = coef + 1.96 * se)) +
  geom_point(position = position_dodge(width = 0.5)) +
  geom_errorbar(position = position_dodge(width = 0.5), width = 0.5) +
  scale_y_continuous(labels = scales::percent_format(accuracy = 1)) +
  labs(x = "Outcome Variable",
       y = "Access likelihood change in percentage points",
       colour = "Social Group",
       shape = "Social Group"
       ) +
  theme_bw() +
  theme(strip.background = element_rect(fill = NA, color = "black"))

ggsave(here("Manuscript","Figures","improvement_access_scst.png"),
       device = "png", width = 6, height = 6*0.618, units = "in")
